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Abstract 

Markov Chain Monte Carlo (MCMC) methods, such as the Metropolis-Hastings (MH) algorithm, 
are widely used for Bayesian inference. One of the most important challenges for any MCMC method 
is speeding up the convergence of the Markov chain, which depends crucially on a suitable choice of 
the proposal density. Adaptive Rejection Metropolis Sampling (ARMS) is a well-known MH scheme 
that generates samples from one-dimensional target densities making use of adaptive piecewise linear 
J> proposals constructed using support points taken from rejected samples. In this work we pinpoint a 

crucial drawback of the adaptive procedure used in ARMS: support points might never be added inside 
regions where the proposal is below the target. When this happens in many regions it leads to a poor 
performance of ARMS, and the sequence of proposals never converges to the target. In order to overcome 
this limitation, we propose two alternative adaptive schemes that guarantee convergence to the target 
distribution. These two new schemes improve the adaptive strategy of ARMS, thus allowing us to simplify 
the construction of the sequence of proposals. Numerical results show that the new algorithms outperform 
the standard ARMS and other techniques. 

Index Terms 

Markov Chain Monte Carlo (MCMC) methods; Metropolis-Hastings (MH) algorithm; Adaptive 
Rejection Metropolis Sampling (ARMS). 
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I. Introduction 

Bayesian inference techniques and their implementation by means of sophisticated Monte Carlo (MC) 
statistical methods, such as Markov chain Monte Carlo (MCMC) and sequential Monte Carlo (SMC) 
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approaches (also known as particle niters), has become a very active area of research over the last years 
IGilks et al.| |1995a| |Gamerman| \TWJ\ |Liu] [20041 |Robert and Casella| [20041 pang et al.[ |20T0[ . Monte 
Carlo techniques are very powerful tools for numerical inference, stochastic optimization and simulation 
of complex systems | |Devroye[ |1986[ [FearnheadJ |1998[ |Doucet et al.[ |2001[ , |Jaeckel[ |2002| , |Robert and 
|CasellaH2004l . 



Rejection sampling (RS) [Devroyej 1986| von Neumann 1951 1 and the Metropolis-Hastings (MH) 



algorithm are two well-known classical Monte Carlo methods for universal sampling [Metropolis and 
Ulaml |1949[ [Metropolis et al.[ |1953[ |Hastings[ |1970| |Liu[ |2004| . Indeed, they can be used to generate 



samples (independent with the RS and correlated with MH) from virtually any target probability density 
function (pdf) by drawing from a simpler proposal pdf. Consequently, these two methods have been 
widely diffused and applied. Unfortunately, in some situations they present certain important limitations 
and drawbacks as we describe below. 

In the RS technique the generated samples are either accepted or rejected by an adequate test of the 
ratio of the target, p(x), and the proposal density, tt(x), with x G V C Rj^An important limitation of 
RS is the need to analytically establish an upper bound M for the ratio of these two pdfs, M > |^ 
or equivalently Mir(x) > p{x), which is not always an easy task. Moreover, even using the best bound, 
i.e. M = sup x p(x) /tt(x), the acceptance rate can be very low if the proposal is not similar to the target 
pdf. On the oher hand, in the MH algorithm, depending on the choice of the proposal pdf, the correlation 
among the samples in the Markov chain can be very high [Liu 2004 Liang et aL] 2010| Martino and 



Mfguez, 2010]. Correlated samples provide less statistical information and the resulting chain can remain 
trapped almost indefinitely in a local mode, meaning that convergence will be very slow in practice. 
Moreover, determining how long the chain needs to be run in order to converge is a difficult task, leading 
to the common practice of establishing a large enough burn-in period to ensure the chain's convergence 



Gilks et al. 


1995a 


Gamerman 


1997 


Liu 2004 



In order to overcome these limitations several extensions of both techniques have been introduced 
[ Devroyej 1986 Liu 2004 Liang et al. 2010], as well as methods combining them [Tierney, 1991, 



1994]. Furthermore, in the literature there is great interest in adaptive MCMC approaches that can speed 
up the convergence of the Markov chain [Haario et al. 20011 Gasemyr) |2003| Andrieu and Moulines) 
20061 lAndrieu and Thorns! [20081 IHolden et all [20091 IGriffin and Walker! |20TT| . Here we focus on two 



'Note that we assume, as is commonly done in the literature, that both pdfs are known up to a normalization constant, and 
we denote as p (x) oc p(x) and tt(x) oc tt(x) the normalized target and proposal pdfs. 
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of the most popular adaptive extensions, pointing out their important limitations, and proposing a novel 
adaptive technique that can overcome them and leads to a much more efficient approach for drawing 
samples from the target pdf. 

One widely known extension of RS is the class of adaptive rejection sampling (ARS) methods |Gilks| 
1992| Gilks and Wild 1992 Gilks et al. 1997 1, which is an improvement of the standard RS technique 



that ensures high acceptance rates with a moderate and bounded computational cost. Indeed, the standard 



ARS algorithm of [Gilks and Wild 1992 1 yields a sequence of proposal functions, {^t{x)}ten, that 
converge towards the target pdf when the procedure is iterated. The basic mechanism of ARS is the 
following. Given a set of support points, <S t = {s\, s mt }, ARS builds a proposal pdf composed of 
truncated exponential pdfs inside the intervals (sj, Sj+i] (0 < j < mt), considering also the open intervals 
(— oo,si] (corresponding to j = 0) and [s mt ,+oo) (associated to j = m t ) when T> = M. Then, when a 
newly proposed sample x' , drawn from 7r t (x), is rejected, it is always added to the set of support points, 
St+i = St U {x'}, which are used to build a refined proposal pdf, m+i(x), for the next iteration. 

Note that the proposal density used in ARS becomes quickly closer to the target pdf and the proportion 
of accepted samples grows. Consequently, since the proposal pdf is only updated when a sample is rejected 
and the probability of discarding a sample decreases quickly to zero, the computational cost is bounded 
and remains moderate, i.e. the computational cost does not diverge due to this smart adaptive strategy 
that improves the proposal pdf only when and where it is needed. Unfortunately, this algorithm can only 
be used with log-concave target densities and hence also unimodal. Therefore, several generalizations 
have been proposed | |Gilks et al| |1995b| |H5rmann] |1995| |Evans and Swartz| |1998| |Goriir and Teh| |2011 



Martino and Migiiez] |2011[ handling specific classes of target distributions or using jointly a MCMC 



approach. 

Indeed, for instance, another popular technique that combines the ARS and MH approaches in an 
attempt to overcome the limitations of both methodologies is the Adaptive Rejection Metropolis Sampling 



(ARMS) [Gil ks et al.[ |1995b[ . ARMS extends ARS to tackle multimodal and non-log-concave densities 
by allowing the proposal to remain below the target in some regions and adding a Metropolis-Hastings 
(MH) step to the algorithm to ensure that the accepted samples are properly distributed according to 
the target pdf. Note that introducing the MH step means that, unlike in the original ARS method, the 
resulting samples are correlated. The ARMS technique uses first an RS test on each generated sample 
x' and, if this sample is initially accepted, applies also an MH step to determine whether it is finally 
accepted or not. If a sample x' is rejected in the RS test (hence, imperatively TVt(x') > p(x'), as otherwise 
samples are always initially accepted by the RS step), then it is incorporated to the set of support points, 
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St+i = St U {x'}, that is used to build a better proposal pdf for the next iteration, Tr t +i(x), exactly as in 
ARS. On the other hand, when a sample is initially accepted, after the MH step the set of support points 
is never updated, i.e. St+i = St, meaning that the proposal pdf is not improved and 7r t+ i(x) = TTt(x). 
This is the crucial point w.r.t. the performance of ARMS. In general, if proper procedures for 



constructing irt(x) given a set of support points, St, such as the ones proposed in |Gilks et al. 1995b 



Meyer et al. [2008 ] are adopted, the proposal will improve throughout the whole domain x G V. However, 



inside the regions of the domain T> where nt(x) < p(x), new support points might never be added. 
Consequently, the convergence of the sequence {iTt(x)}t£n to p(x) cannot be guaranteed, and it may not 
occur in some cases, as shown in Section [TV] through a simple example. Due to this problem and the 
correlation among the samples generated, when the target pdf is multimodal the Markov chain generated 
by the ARMS algorithm tends to get trapped in a single mode despite the (partial) adaptation of the 
proposal pdf (see e.g. example 2 in [Martino and MfguezJ 2010| |). 

It is important to remark that this drawback is caused by a structural problem of ARMS: the adaptive 
mechanism proposed for updating the proposal pdf is not complete, since it never adds support points 
inside regions where irt(x) < p(x). Therefore, even though a suitable procedure for constructing the 
proposals (see e.g. | Gilks et ah] 1 1995b) |Meyer et al. 2008 1) can allow them to change inside these regions 
(eventually obtaining nt(x) > p{x) in some cases), the convergence cannot be guaranteed regardless of 
the procedure used to build the proposal densities. For this reason, the performance of the standard ARMS 
depends critically on the choice of a very good way to construct the proposal pdf, 7Tt(x), that attains 
Kt{x') > p(x') almost everywhere, implying that ARMS tends to be reduced to the ARS algorithm. 
Indeed, note that if the procedure used to build the proposal produces 7ro(x) < p(x) Vx G V, then the 
proposal is never improved, i.e. nt(x) = ttq(x) for all t G N, resulting in no adaptation at all. 

This structural problem of the ARMS approach can be solved in a trivial way: adding a new support 
point each time that the MH step is applied. Unfortunately, in this case the computational cost of the 
algorithm increases rapidly as the number of support points diverges. Moreover, a second major problem 
is that the convergence of the Markov chain to the invariant density cannot be ensured, as partially 
discussed in jGilks et al] |1995b] and in pGilks et al!] |1997Q for the case of ARMS-within-Gibbs (see 
also [Liang et al. 2010| Chapter 8] for further considerations). 

In this work we present two enhancements of ARMS that guarantee the convergence of the sequence 
of proposal densities to the target pdf with a bounded computational cost, as in the standard ARS and 
ARMS approaches, since the probability of incorporating a new support point quickly decreases to zero. 
Moreover, the two novel techniques fulfill one of the required condition to assure the convergence of 
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the chain, the so-called diminishing adaptation (see e.g. I Liang et al. 2010[ Chapter 8], [Haario et al. 



2001 1 Roberts and Rosenthal} 2007| ). The first one is a direct modification of the ARMS procedure that 



allows us to incorporate support points inside regions where the proposal is below the target with a 
decreasing probability as the proposal converges to the target. The second one is an adaptive independent 
MH algorithm that learns from all past generated samples except for the current state of the chain, thus 
also guaranteeing the convergence of the chain to the invariant density, as shown in [ Gasemyr] 2003 



Holden et al. 2009 1. Furthermore, these new strategies also allow us to reduce the complexity in the 



construction of the proposal pdfs (thus reducing both the effort of writing the code and the computational 
cost of the resulting algorithm), since they do not require that TTt(x) > p{x) almost everywhere as in the 
standard ARMS algorithm in order to guarantee the smooth running of the adaptive mechanism and hence 
to improve the proposal pdf everywhere in the whole domain V. We exemplify this point introducing 
simpler procedures to construct the sequence of proposal densities and illustrating their good performance 
through numerical examples. 

We remark that, compared to other adaptive Metropolis-Hastings techniques available in literature 
| Warnes[ |200T| |Haario et al] [200T] [20061 |Cai et al.| [2008] pang et al.| |20"T0| , the two new schemes 
provide a better performance. For instance, two alternative adaptive schemes used to draw samples from 
one-dimensional target pdfs are the adaptive triangle Metropolis sampling (ATRIMS) and the adaptive 



trapezoid Metropolis sampling (ATRAMS) proposed in pCai et aL] |2008| |. However, even though the 
proposal pdfs are adaptively improved in both of them (ATRIMS and ATRAMS), the sequence of 
proposals does not converge to the target density, unlike the adaptive strategies proposed in this work. 
Regarding other adaptive MH algorithms | |Warnes| |20"0lj |Haario et al"j [200T] [20061 |J- M. Keith~and 



Sofronov 2008 [ |Hanson et al.} 201 1[ only certain parameters of the proposal pdf (which follows a 



parametric model with a fixed analytical form) are adjusted and optimized, whereas here we improve 
the entire shape of the proposal density, which becomes closer and closer to the shape of the target 
density. Finally, another advantage of our schemes is that the proposed algorithms eventually become 
standard ARS techniques when the constructed proposal pdf is always above the target, i.e. when 
T^t{x) > p(x) Vx 6 V, thus producing independent samples with acceptance rate that quickly becomes 
close to one. For all these reasons, drawing from univariate pdfs, our techniques obtain better performance 
than the other techniques available in literature at the expense of a moderate increase in the computational 
cost. 

The rest of the paper is organized as follows. Background material is presented in Section [IT] Then we 
review and discuss certain limitations of ARMS in Sections III and IV In Section [V] we present the two 
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novel techniques, whereas alternative procedures to construct the proposal pdf are described in Section 



VI Finally, numerical simulations are shown in Section VII and conclusions are drawn in Section VIII 



II. Background 

A. Notation 

We indicate random variables with upper-case letters, e.g. X, while we use lower-case letters to denote 
the corresponding realizations, e.g. x. Sets are denoted with calligraphic upper-case letters, e.g. 1Z. The 
domain of the variable of interest, x, is denoted as V C R. The normalized target PDF is indicated 
as p (x), whereas p{x) = c p p (x) represents the unnormalized target. The normalized proposal PDF is 
denoted as tt(x), whereas tt(x) = c 7r fr(x) is the unnormalized proposal. For simplicity, we also refer 
to the unnormalized functions p(x) and n(x), and in general to all unnormalized but proper PDFs, as 
densities. 

B. Adaptive rejection sampling 

The standard adaptive rejection sampling (ARS) algorithm | Gilks| 1992| [Gilks and Wild 1992 Gilks 



et al. 1997| enables the construction of a sequence of proposal densities, {^t(x)} ten , tailored to the 
target density, p (x) oc p{x). Its most appealing feature is that each time that a sample is drawn from a 
proposal, TTt(x), and it is rejected, this sample can be used to build an improved proposal, Tr t +i(x), with a 
higher mean acceptance rate. Unfortunately, the ARS method can only be applied with target pdfs which 
are log-concave (and thus unimodal), which is a very stringent constraint for many practical applications 



| Gilks et al. 1995b Martino and Mi'guez 2011 1. In this section we briefly review the ARS algorithm, 
which is the basis for the ARMS method and the subsequent improvements proposed in this paper. 

Let us assume that we want to draw samples from a target pdf, p {x) oc p(x), with support PCI, 
known up a normalization constant. The ARS procedure can be applied when p{x) is log-concave, i.e. 
when 

V{x) 4 \og(p{x)), (1) 

is strictly concave \/x € V C R. In this case, let 

St = {si,s 2 , ■ ■ ■ ,s m J C V (2) 

be the set of support points at time t, sorted in ascending order (i.e. si < . . . < s mt ). Note that the 
number of points m t can grow with the iteration index t. From St a piecewise-linear function Wt(x) is 
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constructed such that 

W t (x) > V(x), (3) 
for all x G V and Vi G N. Different approaches are possible for constructing this function. For instance, 



if we denote by Wk{x) the linear function tangent to V(x) at Sk G St | Gilks and Wild 1992 1, then a 
suitable piecewise linear function Wt(x) is: 



Wt(x) = mm{w\(x), . . . ,w mt (x)} for i£P. 



(4) 



In this case, clearly we have Wi(x) > F(x), Vx G V and Vi G N, as shown in Figure [jj a). However, 
it is also possible to construct Wt{x) without using the first derivative of V(x), which is involved in 
the calculation of the tangent lines used in the previous approach, by using secant lines [Gilks 1992| |. 
Indeed, defining the intervals 

Z = (-oo,si], Xi = (si,s 2 ], ■■■ - (sj, s j+i], Zm t - ( s m t -i,s mt ], T-mt - (s mt ,+oo) (5) 

and denoting as Li^ + \{x) the straight line passing through the points (sj, V{si)) and (sj+i, V(si + i)) for 
1 < i < m t , a suitable function Wt{x) can be expressed as 

L 1)2 (x), x G X = (-00, si]; 

^2,3^), iGl 1 = (si,s 2 ]; 

W t (x) = <[ min^i^x), L i+1)i+2 (x)}, x Elj = ( Sj , s j+1 ], 2 < j < m t - 2; (6) 

^'nii - 2, rot -• 1 (jE) ; *E G X mt _x — (s mf _i,S mt ], 

L'm t — l,m t (^)j 2- G 2m t — (^m t i "^Oo). 

It is apparent also in this case that Wt(x) > V(x), \/x £ V and Vi G N, as shown in Figure [jjb). 
Therefore, building Wt(x) such that Wt(x) > V(x), we have in both cases 

7T t (x) = exp(V^(x)) > p(x) = exp{V(x)}. (7) 

Since Wt(x) is a piecewise linear function, we obtain an exponential-type proposal density, 7ft (x) oc -Kt (x). 
Note also that knowledge of the area below each exponential piece, ki = J xeX , ex.p(Wt(x))dx for 
< i < m t , is required to generate samples from 7f 4 (x), and, given the functional form of 7r 4 (x), it 
is straightforward to compute each of these terms, as well as the normalization constant, l/c n with 
c tt = J x£V n(x)dx, since = k + k 2 + ... + k n . 

Therefore, we conclude that, since ir t (x) is a piecewise-exponential function, it is very easy to 
draw samples from 7ft(x) = — TTt(x) and, since 7T((x) > p(x), we can apply the rejection sampling 
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principle. Moreover, when a sample x' drawn from 7ft (x) is rejected, we can incorporate it into the 
set of support points, i.e. St+i = St U {x'} and m t +i = mj + 1. Then, we can compute a refined 
approximation, Wt+i(x), following one of the two approaches described above, and a new proposal 
density, 7r t +i(x) = exp(W/ 7 t+ i(x)), can be constructed. Table [i] summarizes the procedure followed by 
the ARS algorithm to generate N samples from a target p(x). 



TABLE I 

Adaptive Rejection Sampling Algorithm. 



1. Start with i = 0, t = 0, mo = 2, <So = {si, S2} where si < S2 and [si, S2] contains the mode. 
Let N be the number of desired samples from p (x) oc p(x), 

2. Build the piecewise-linear function Wt(x) as in Eq. i4i or 161, and as shown in Figure jl|a) orjljb). 

3. Sample x' from nt(x) oc ir t (x) = exp(Wt(x)), and it' from W([0, 1]). 

4. If u < then accept — x' and set St+i = 5t (mt+i = mt), i = i + 1. 



5. Otherwise, if it' > , A , then reject x', set 5t+i = St U {a;'} (mt+i = m t + 1). 

6. Sort 5t+i in ascending order and set t = t + 1. If j = A, then stop. Else, go back to step 2. 



It is important to observe that the adaptive strategy followed by ARS includes new support points 
only when and where they are needed, i.e. when a sample is rejected, indicating that the discrepancy 
between the shape of the target and the proposal pdfs is potentially large. This discrepancy is measured 
by computing the ratio p(x)/irt(x). Since nt(x) approaches p{x) when t — > 00, the ratio p(x)/nt(x) 
becomes closer and closer to one, and the probability of adding a new support point becomes closer 
and closer to zero, thus ensuring that the computational cost remains bounded. More details about the 
convergence of ARS are provided below. 

C. Convergence 

Note that every time a sample x' drawn from nt(x) is rejected, x' is incorporated as a support point in 
the new set St+i = St U {x'}. As a consequence, a refined lower hull Wt+i{x) is constructed, yielding 
a better approximation of the system's potential function, V{x). In this way, -K t +i{x) = exp(Wf+i(x)) 
becomes "closer" to p(x) and the mean acceptance rate can be expected to become higher. More precisely, 
the probability of accepting a sample x G V drawn from irt(x) is 

a t (x) 4 P^l (8) 



9 




Fig. 1. Two possible procedures to build a suitable piecewise linear function Wt(x) > V(x) when V(x) is concave: (a) using 
the first derivative of V(x) (i.e. tangent lines) at the support points (with 3 support points, for instance); (b) using secant lines 
(with 5 support points, for instance). 



and we define the acceptance rate at the t-th iteration of the ARS algorithm, denoted as dt, as the expected 
value of at(x) with respect to the normalized pdf nt(x), i.e., 

at = E[at{x)} = f a t (x)TTt(x)dx = f ^ \ iit{x)dx = — f p(x)dx = — < 1, (9) 



v Jv 



where l/c n and l/c p are the proportionality constants for irt(x) and p(x) respectively, i.e. 

c-k = / irt(x)dx, 
Jv 

and 



c p = p(x)dx, 
Jv 

and we remark that ^ < 1 because TTt(x) > p(x) Vx G V. Hence, from (|9]) we conclude that 
dt = 1 4=> c n = c p . Equivalently, defining the distance between two curves D^\ p (t) as 



D^\p{t)= / \7Tt(x)-p(x)\dx= / |exp(W t (x))-exp(^(x))|dx, (10) 
Jv Jv 



then we have dt = 1 4=> D n \ p (t) = 0. Note that this distance D n \ p (t) measures the discrepancy between 
the proposal irt(x) and the target p{x) pdfs. In particular, if D^ p (t) decreases the acceptance rate dt = — 
increases, and, since exp(Wt(x)) > exp(V(x)) \/x G V, D n \ p (t) = if, and only if, Wt(x) = V(x) 
almost everywhere. Equivalently, d t = 1 if, and only if, TT t (x) = p(x) almost everywhere. 
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III. Adaptive Rejection Metropolis Sampling Algorithm 

The disadvantage of the classical ARS technique is that it can only draw samples from univariate log- 
concave densities. In general, for non-log-concave and multimodal pdfs it is difficult to build a function 



Wt(x) that satisfies the condition Wt{x) > V(x) for all x G T>. Recognizing this problem, in [Gilks 



et al.[ |1995b[ an extension of ARS is suggested to deal with pdfs that are not log-concave by appending a 



Metropolis-Hastings (MH) algorithm step [Hastings 1970| Metropolis et al. 1953| Metropolis and Ulam 



1949]. The adaptive rejection Metropolis sampling (ARMS) first performs an RS test, and the discarded 



samples are used to improve the proposal pdf, as in the standard ARS technique. However, if the sample 
is accepted in the RS test, then the ARMS adds another statistical control using the MH acceptance 
rule. This guarantees that the accepted samples are distributed according to the target pdf, p(x), even if 
Wt(x) < V{x). Therefore, unlike ARS, the algorithm produces a Markov chain and the resulting samples 
are correlated. Finally, we remark that the ARMS technique can be seen as an adaptive generalization of 
the rejection sampling chain proposed in [Tierney, 1991[ 1994 1 . 

The ARMS algorithm is described in detail in Table [n] The time index k denotes the iteration of the 
chain whereas t is the index corresponding to evolution of the sequence of proposal pdfs {Ttt(%)}teN- 

The key stages are steps 4 and 5. On the one hand, in step 4 of Table [nj when a sample is rejected 
by the RS test (it is possible if and only if irt(x') > p(x')), ARMS adds this new point to the set St 
and uses it to improve the construction of Wt+i(x) and go back to step 2. On the other hand, when a 
sample is initially accepted by the RS test (clearly this always happens if TTt(x') < p{x')), ARMS enters 
in step 5 and uses the MH acceptance rule to determine whether it is finally accepted or not. However, 
notice that ARMS never incorporates this new point to St, even if it is finally rejected by the MH step. 
Observe also that, if nt(x) > p(x) Mx £ V and Vt G N, then the ARMS becomes the standard ARS 
scheme, since the proposed sample is always accepted in step 5. 

Moreover, since 7Tt(x) depends on the support points, a more rigorous notation would be ir t (x\St). 
However, since St never contains the current state of the chain, then the ARMS can be considered an 
adaptive independent MH algorithm and, for this reason, the dependence on St is usually omitted from 
the notation. Indeed, the ARMS can be also seen as a particular case of the auxiliary variables method 
in [Besag and Green| |1993[ . 



Finally, we observe that, although in [Gilks et al. 1995b | the ARMS algorithm seems to be tightly 
linked to a particular mechanism to construct the sequence of the proposal pdfs, in fact these two issues 
can be studied separately. Therefore, in the following section we describe the two procedures proposed 
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in literature to construct the function Wt(x) in a suitable way [Gilks et al. 1995b Meyer et al. 2008]. 



TABLE II 

Adaptive Rejection Metropolis Sampling Algorithm. 



1. Start with k — (iteration of the chain), t — 0, So = {si, . . . , s mo }, with si < S2 < ■ ■ ■ < s mo , 
and choose a value x$. Let N be the required number of iterations of the Markov chain. 

2. Build a function Wt(x) using the set of points St and one of the two procedures described in Section 



Gilks et al. 


1995b 


Meyer et al. 


2008 



III-A 



3. Sample x' from nt(x) oc n t (x) — cxp(Wt(x)), and u from W([0, 1]). 

4. If u > J^ry, then discard x', set St+i — St U {x'}, mt+i — m t + 1, sort St+i in ascending order, 
update t = t + 1, and go back to step 2. 



5. Otherwise, i.e. if u' < 

a 



Tr t (x') 

min 



jt, set xt+i = x' with probability 



N p(z') min[p(3; fc ),7rt(3: fc )] "| 

or set Xk+i ~ Xk with probability 1 — a. Moreover, set St+i ~ St, rrit+i = m t , t - 
6. If k = N — 1, then stop. Else, go back to step 2. 



£ + landfc = fc + l. 



A. Procedure to build Wt(x) proposed in [Gilks et al. 1995b^ 

Consider again a set of support points St = {si, ...,s mt } and the m t + 1 intervals 2o = (— oo,si], 
Ij = (sj,Sj + i\, for j = l,...,mt — 1 and X mt = (s mt ,+oo). Moreover, let us denote as Ljj + i(x) the 
straight line passing through the points (sj,V(sj)) and (sj+i,V(sj+i)) for j = l,...,m t — 1, and also 
set 

L_i i0 (aj) = L 0) i(x) = Li )2 (x), 

Lrn t ,m t +l( x ) = ^m t +l,m t +2 {x) = L mt —l >rn , t (x). 



In the standard ARMS procedure introduced in | Gilks et al.| 1995b |, the function Wt{x) is piecewise 
linear and defined as 

W t (x) = max [Li ] i +1 (x),mm[Li^i > i(x),L i+ i ] i + 2(x)]] for Xj = (sj, s i+1 ], (11) 
with i = 0, . . . , 777.4. The function Wt(x) can be also rewritten in a more explicit form as 

f 

L 1<2 (x), x €Iq = (-cx5,sj; 

max{£i )2 (a;),L 2i 3(x)} , ieli = 0i,s 2 ]; 

Wt(a;) = { max{Ljj + i(x),min{Lj-ij(x),Lj + i jj+ 2(x)}}, x elj = (sj,s j+1 ], 2 < j < m t - 2; (12) 

max{L mt _i ;mt (x), L mt _ 2 ,m t -i{x)} , x el mt -i = (s mt -i,s mt ]\ 



L 



mt-l 



x ^ X mt (s mt , -(-oo). 
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Figure [2] illustrates this construction for a generic log-pdf V(x) with (a) 5 support points and (b) 6 
support points. 




(a) 



(b) 



Fig. 2. Example of piecewise linear function Wt(x) built using the procedure described in | Gilks et al. 1995b) : (a) with 5 
support points and (b) with 6 support points. Note that V(x) is non-concave and Wt{x) can fall below V(x). 



It is important to remark that, with this construction, the addition of a new point inside an interval 
Ti = (sj,Sj + i] also affects the adjacent regions, Zj_i = (sj_i,Sj] and Zj + i = (si+i , Si+2] • Therefore, if 
the initial proposal pdf iro(x) fulfills the inequality ttq(x) > p(x) in some subset C C V, the proposal 
pdf can theoretically be improved in the entire domain V. However, in practice the convergence of the 
proposal to the target pdf can be very slow in some cases, and moreover it cannot be ensured in general 



(see Section [TV] for further considerations). 

Note also that, if the target pdf p(x) is log-concave, then 

max[Ljj + i(x),min [Lj-ij(x), L j+ ij +2 (x)]] = mm[Lj-ij(x), L j+ ij +2 (x)] , 

for < j < m t , hence this procedure is reduced to Eq. ([6]) of the standard ARS technique. It is important 



to observe that, although Wt(x) can be expressed in the compact form of Eq. (Ill, drawing samples from 
the proposal, TTt(x) oc exp(Wt(x)), requires knowing the equation of each straight line that composes 
Wt(x), and also the corresponding definition interval (see Figure [2]). Hence, the complexity of drawing 



from TTt(x) is much larger than what equations ( 11 ) or (12) might suggest, since it entails the computation 
of all the intersection points between all the contiguous straight lines (see Figure [2]). 

Simpler constructions without this requirement could be proposed (see, for instance, Section |Vl]), but 
the efficiency of the standard ARMS scheme would be severely compromised, since the good performance 
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of ARMS depends critically on the construction of the proposal, as discussed in the following section 
and shown in the simulations (see Section |VII| ). 



Another possible and more complex procedure for building Wt(x) was proposed in [Meyer et al. 



2008|, involving quadratic approximations of the log-pdf V(x) when it is possible, so that the proposal 
is also formed by truncated Gaussian pdfs. Indeed, parabolic pieces passing through 3 support points 
are used if only if the resulting quadratic function is concave (in order to obtain truncated Gaussian as 
proposal pdf in the corresponding interval). Otherwise, linear functions are used to build Wt(x). The main 



advantage of this alternative procedure [Meyer et al. 2008 1 is that it provides a better approximation of 



the function V(x) = \og(p(x)). However, it does not overcome the critical structural limitation of the 
standard ARMS technique that we explain in the sequel. 

IV. Structural limitation in the adaptive procedure of ARMS 

Although ARMS can be a very effective technique for sampling from univariate non-log-concave pdfs, 
its performance depends critically on the following two issues: 

a) The construction of Wt(x) should be such that the condition Wt(x) > V(x) is satisfied for 
most intervals x G Ij as possible, with j = 0, 1, . . . , m t . That is, the proposal function 
TTt(x) = exp(— Wt(x)) must stay above the target pdf p(x), inside as many intervals as possible, 
covering as much of the domain V as possible. In this case the adaptive procedure works almost 
in the entire domain V and the proposal density can be improved virtually everywhere. Indeed, if 
iTt(x) > p(x) for all x G V and for all t G N, the ARMS is reduced the classical ARS algorithm. 

b) The addition of a support point within an interval, Ij = (sj, Sj + i], with j G {0, m t }, must entail 
an improvement of the proposal pdf inside other neighboring intervals when building Wt+\(x). This 
allows the proposal pdf to improve even inside regions where irt(x) < p(x) and a support point 
can never be added at the t-th iteration, since we could have 7r t+ i(x) > p(x) by adding a support 
point in an adjacent interval. For instance, in the procedure described in Section [Til- A [Gilks et al. , 
1995b |, when a support point is added inside Ij, the proposal pdf also changes in the intervals 



Ij-i and Ij+i. Consequently, the drawback of not adding support points within the intervals where 
TTt(x) < p(x) is reduced, but may not completely eliminated, as shown through an example below. 
In any case, it is important to remark that the convergence of the proposal irt(x) to the target pdf 
p(x), cannot be guaranteed using any suitable construction of Wt(x) except for the special case where 
Wt(x) > V{x) \/x G V and Vt G N, where ARMS is reduced to ARS. This is due to a fundamental 
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structural limitation of the ARMS technique caused by the impossibility of adding support points inside 
regions where ir t (x) < p{x). 

Indeed, it is possible that inside some region C C V where TTt(x) < p(x), we obtain a sequence of 
proposals {ir t+1 (x), 7r i+2 (x), . . . , 7r t+r (x)} such that 7r m (x) < p(x), TT t+2 (x) < p(x), ir t+T (x) < 
p{x) for an arbitrarily large value of r, and the discrepancy between the proposal and the target pdfs 
cannot be reduced after an iteration t = r, clearly implying that the proposal does not converge to the 
target for 

In order to illustrate this structural problem of ARMS we provide below an example of an even more 
critical situation where ~k%{x) = ir t +i(x) = 7r t + 2 (x) • • • for all x £ Ij, i.e. the proposal pdf does not 
change within an interval Ij C V. Consider a multi-modal target density, p {x) oc p{x) = exp(V(x)), 
with function V(x) as shown in Figure |3ja). We build Wt(x) using 5 support points and the procedure of 



Section |iIT"A"| jGilks et"aT}|1995b[ . Note that we have W t (x) < V(x) for all x in the interval I 2 = (s 2 , S3], 
as shown in Figure [3^a), where the dashed line depicts the tangent line to V{x) at S3. Figures |5Jb) and 
J3jc) show that we can incorporate new support points (S4 in Figure [3jb) and s 2 in Figure [3jc)) arbitrarily 
close to the interval Z3 (I 2 in Figures |5|a)-(b)) without altering the construction of Wt{x) within Z3 (X 2 
in Figures [3ja)-(b)). The reason is described in the sequel. 

Consider Figure a). The construction of Wt(x) for all x G X 2 depends on the straight lines L\^{x) 
passing through {si,V(s\)) and (s 2 ,F(s 2 )), L 2 ^{x) passing through (s 2 ,F(s 2 )) and (s 3 , V(s 3 )), and 
L3,4(x) passing through (s3,y(s3)) and (S4, V(s4)). Hence, Eq. ( [TT] ) within this interval can be written 
as 

W t (x) = max (L 2j3 (a;), min {L lt2 (x), L 3A (x)} }, for x € 1 2 = (s 2 , s 3 ]. (13) 

Looking back at Figure [5|a) it becomes apparent that min {Li 2 (x), L3 4(x)} = £3,4(3;) and 
max{L 2: 3(x), L^^{x)} = L 2 ^(x). Therefore, Wt{x) = L 2 ^(x) for all x G 1 2 = [s 2 ,S3], and this 
situation does not change when new support points are added inside the contiguous intervals, as shown 
in Figures |3jc)-(d). Indeed, consider now the limit case where two points are incorporated arbitrarily close 
to s 2 and S3, so that we can consider the new support points located exactly at s 2 and S3. In this extreme 
situation the secant lines of the adjacent intervals become tangent lines, as shown in Figure [3jd), and 
the minimum between the two tangent lines is represented by the straight line tangent to S3. Moreover, 
note that this tangent line stays always below the secant line L 2j 3(x) passing through (s 2 ,V(s 2 )) and 
(S3, V(s3)), meaning that Wt(x) = L 2 ^(x) even in this extreme case. 
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(a) 



(b) 





(c) 



(d) 



Fig. 3. Example of a critical structural limitation in the adaptive procedure of the ARMS, (a) Construction of Wt{x) with 5 
support points. Within I2 = (S2, S3] we have Wt(x) < V(x). (b)-(c) Adding new support points inside the contiguous intervals 
the construction of Wt(x) does not vary within X2 (I3 in Figure (c)). (d) The secant line L2,i(x) passing through (S2, V(s2)) 
and (S3, V(s3)), and the two tangent lines to V(x) at S2 and S3. 



Therefore, with the adaptive procedure used by the standard ARMS technique to build the proposal 
pdf, the function Wt(x) could remain unchanged in a subset of the entire domain V that contains a 
mode of the target pdf, as shown in the previous example. Hence, the discrepancy between the proposal 
iTt(x) = exp(Wt(x)) and the target p(x) = exp(V(x)), is not reduced at all, as t — > +00. 

A trivial solution for this drawback of the standard ARMS algorithm could be adding new support 
points within the set St each time that the MH step in the ARMS (step 5) is applied. Unfortunately, from 
a practical point of view the first problem of this approach is the unbounded increase in computational 
cost since the number of points in St grows indefinitely. Another major problem from a theoretical point 
of view is that the convergence of the Markov chain to the invariant pdf cannot be ensured in this case, 
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since the current state of the chain could be contained in St (see [Liang et al. 2010[ Chapter 8], |Holden 



et al.| |2009| | for a discussion of this issue). 

For these two reasons in the following section we propose two variants of the standard ARMS technique 
that ensure the convergence of the Markov chain to the target density while keeping the computational 
cost bounded and low. 

V. Variants of the ARMS algorithm 

In this section we describe two alternative strategies to improve the standard ARMS algorithm. We 
denote these two variants as A 2 RMS and IA 2 RMS where the A 2 emphasizes that we incorporate an 
additional adaptive step to improve the proposal density w.r.t. the standard ARMS. Note that in this 
section we use the more rigorous notation for the proposal, 7Tt(x\St) instead of simply vr t (x), for clarity 
with respect to the convergence of the Markov chain. 

A. First scheme: A 2 RMS 

A first possible enhancement of the ARMS technique is summarized below. The basic underlying 
idea is enhancing the performance of ARMS by introducing the possibility of adding a sample to the 
support set even when it is initially accepted by the RS test. The procedure followed allows A 2 RMS 
to incorporate support points inside regions where nt(x\St) < p(x) in a controlled way (i.e. with a 
decreasing probability as the Markov chain evolves), thus ensuring the convergence to the target and 
preventing the computational cost from growing indefinitely. 

The steps taken by the A 2 RMS algorithm are the following: 

1. Set k = (iteration of the chain), t = 0, choose an initial value xq and the time to stop the 
adaptation, K. Let N > K be the needed number of iterations of the Markov chain. 

2. Given a set of support points, St = {s±, . . . , s mt } such that s± < S2 < ■ ■ ■ < s mt , build an 
approximation Wt(x) of the potential function V(x) = log(p(x)) using a convenient procedure 
(e.g. the ones described in [Gilks et ah] 1995b[ Meyer et al.[ |2008] or the simpler ones proposed in 



Section |VJ. 

3. Draw a sample x' from n t (x\St) oc 7Tt(x\St) = exp(Wi) and another sample from u' ~ W([0, 1]). 

4. If v! > ^v|,s t ) > then discard x' , set St+i = St U {x'}, mt+i = writ + 1, sort St+i in ascending 
order, update t = t + 1, and go back to step 2. 

5. If u' < ffiL , then: 

— 1T t (x'\St)' 
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5.1. Set Xk+i = x' 



■' with probability 



a = min 1 



p(x') mm[p(x k ),7r t (x k \St)} 
' p(x k )mm[p(x'),Tr t (x'\S t )} 



(14) 



or ^fc+i = x k with probability 1 — a. 
5.2. If k < K, draw u 2 ~ W([0, 1]) and if 



«2 > 



set <St+i = StU{x'}, mt+i = rrit + l and sort St+i in ascending order. Otherwise, set St+i = St, 
m t+1 = m t . 
5.3. Update t = t + 1 and k = k + 1. 
6. If fc < N, go back to step 2. 
We remark again that the key point in the A 2 RMS algorithm is that, due to the introduction of step 
5.2 w.r.t. the ARMS, new support points can also added inside the regions of the domain V where 
irt(x\St) < p(x). Note that, when -K t {x'\St) > p(x') and we accept the proposed sample x' in the 
RS test, then we also apply the MH acceptance rule and always accept x' , as a = 1, but we never 
incorporate x' to the set of support points, since ^fc^ > 1 and u 2 ~ U([0,1\). Therefore, if the 
condition 7r t (x\St) > p(x) is satisfied in all the domain x G V, then the A 2 RMS is reduced to the 
standard ARS, just like it happens for ARMS. 

It is important to remark that the effort to code and implement the algorithm is virtually unchanged. 
Moreover, the ratio ^feif^ = p( i/) — hi the step 5.2 is already calculated in the RS control (steps 4 



and 5), hence it is not necessary to evaluate the proposal and the target pdfs again. 

Since it is difficult to guarantee that the Markov chain converge to the target pdf, p {x) oc p(x), we 
have introduced a time K to stop the second adaptation step 5.2. Therefore, theoretically the first K 
samples produced by the algorithm should be discarded. Observe that for k > K the A 2 RMS coincides 
to the standard ARMS. 

The issue with the convergence of the chain is due to the fact that we are incorporating the current state 
xt, into the set of support points St, on which the proposal depends. Therefore, the balance condition 
using the acceptance function in Eq. ( [T4| ) could be not satisfied for k < K and those first K samples could 
be seen just as auxiliary random variables obtained as part of the process to construct a good proposal 
pdf, and thus should be removed from the set of final returned samples. However, it is important to notice 
the following two facts: 

a) It is a common practice with MCMC techniques to remove a certain amount of the firstly generated 
samples in order to diminish the effect of the so called burn-in period. 
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b) We have found out empirically that, even if we set K = N and use all the samples produced by 
the Markov chain, the A 2 RMS algorithm also provides very good results, as we show in Section 



|VII| Indeed, the probability of incorporating new support points quickly tends to zero as t increases, 
effectively vanishing as t — > +00. This is due to the fact that, each time that a new point is added, 
the proposal density is improved and becomes closer and closer to the target pdf in the whole 
domain V. Therefore, the two ratios n j^ St ^ and approach one and the probability of adding 

a new support point becomes arbitrarily close to zero. This implies that the computational cost of 
the algorithm remains bounded, i.e. a very good approximation of the target can be obtained with 
a finite number of support points adaptively chosen. 
Hence, A 2 RMS satisfies the first condition needed to ensure the convergence of the Markov chain to the 
target pdf, known as diminishing adaptation (see JLiang et~aE] |2010| Chapter 8], [Haario et al.| |2001 



Roberts and Rosenthal] |2007| ]). Unfortunately, it is difficult to assert whether the A 2 RMS with K = N 
also fulfills the second condition needed to guarantee the convergence of the chain, called bounded 
convergence | Liang et al. 2010[ Chapter 8]. For this reason, although the good results obtained in the 



numerical simulations with K = N (see Section VII ) lead us to believe that convergence occurs, it may 



be safer to set K < N in order to avoid convergence problems in practical applications. 

Furthermore, in the next section we introduce a second A 2 RMS scheme for which convergence to the 
target can be ensured, since it is an adaptive independent MH algorithm | Gasemyr[ 2003 Holden et al. 
20091 . 



B. Second scheme: Independent A 2 RMS (IA 2 RMS) 

A second possible improvement of ARMS is an adaptive independent MH algorithm that we indicate as 
IA 2 RMS. The main modification of IA 2 RMS w.r.t. A 2 RMS is building the proposal pdf using possibly 
all the generated past samples but without taking into account the current state x% of the chain. The 
IA 2 RMS is described in the following. 

1. Set k = (iteration of the chain), t = and choose an initial value xq. 

2. Given a set of support points, St = {si, • • • , s mt }, such that s\ < S2 < ■ ■ ■ < s mt , build an 
approximation Wt(x) of the function V(x) = log (p(x)) using a convenient procedure (e.g. the ones 
described in [Gi lks etliE] |1995b[ |Meyer et al.| |2008| | or the simpler ones proposed in Section |Vl] ). 



3. Draw a sample x' from nt{x\St) oc irt(x\St) = exp(Wi) and another sample from v! ~ U{[0, 1]). 

4. If u' > - > , then discard x', set 5 t+ i = St U {x'} and mt+i = m t + 1, sort St+i in ascending 
order, update t = t + 1 and go back to step 2. 
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5. If u> < then: 

5.1 Set Xk+i = x' and y = Xk with probability 

. [ p(x')mm\p(x k ),n t (x k \S t )}~\ 

or Xk+i = Xk and y = x' with probability 1 — a. 

5.2 If y is not already contained in St, draw ii2 ~ W([0, 1]), and if 

7rt(j/|«S f ) 

p{y) 

set <St+i = 5tU{y}, m t +i = m* + l and sort <St+i in ascending order. Otherwise, set St+i = St 
and mt+i = mt- 

5.3 Update t = t + 1 and k = k + 1. 

6. If fc < N go back to step 2. 

Note that the main difference of IA 2 RMS w.r.t. A 2 RMS lies in steps 5.1 and 5.2, where an auxiliary 
variable y is introduced to determine whether a new point is added to the set of support points or not 
instead of using directly x' . 

This construction leads to a proposal, 7r t (x\St), which is independent of the current state x t of the 
chain, and in |Gasemyr[|2003[|Holden et al.[|2009| l, |Liang et al.l|2010[ Chapter 8] it is proved that for this 



kind of adaptive independent MH techniques the convergence of the Markov chain to the invariant density 
is ensured. Once more, if the procedure followed to build the proposal pdf results in TT t (x\St) > p(x) 
for all x G V and t € N, then IA 2 RMS is reduced to the standard ARS algorithm, just like A 2 RMS and 
ARMS. 

Moreover, note that also in the IA 2 RMS the two components ir t (y\St) and p(y), in the ratio » 



are already calculated either if y = x' or y = xt for instance in Eq. (15 1, hence it is not necessary to 
evaluate the proposal and the target densities again. 

Finally, note that new statistical information is used in the two improved ARMS algorithms proposed 
(A 2 RMS and IA 2 RMS) to update the proposal density also in the regions of the space where 
7r t (x\St) < p{x). Consequently, these two algorithms are no longer subject to the restrictions of ARMS 
for attaining a good performance, i.e. that Wt{x) > V(x) in most intervals and the improvement of the 
proposal in adjacent regions when a new support point is added. Therefore, as shown in the following 
section, simpler procedures can be used to build the function Wt(x), reducing the overall computational 
cost of the resulting algorithms and making it easier to code them. 



20 



VI. Alternative procedures to build n t (x) 

As discussed before, the improvements proposed in the structure of the standard ARMS allow us to 
use simpler procedures to construct the function Wt(x). 

For instance, a first approach inspired by the one used in ARMS is defining Wt(x) inside the z-th 
interval as the straight line going through (sj, V(si)) and (sj+i, V(si+i)), Li t i+i(x), for 1 < i < mt — 1, 
and extending the straight lines corresponding to X\ and I mt -i towards minus and plus infinity for the 
first and last intervals, Zq and X mt respectively. Mathematically, this can be expressed as 



W t (x) 



Li,2(x), 



L; 



x G Zo = (-oo, si]; 

x eli = (si, Sj+i], 1 < i < mt - 1; 



(16) 



L 



x G I, 



(s mt! +oo). 



This is illustrated in Figure[4ja). Note that, although this procedure looks similar to the one used in ARMS, 
as described by Eqs. (|TTj» and (12 1, it is much simpler in fact, since there is not any minimization or 
maximization involved, and thus it does not require the calculation of intersection points to determine 
when one straight line is above the other. 

Moreover, an even simpler procedure to construct Wt(x) can be devised from ( fT6] ): using a piecewise 
constant approximation with two straight lines inside the first and last intervals. Mathematically, it can 
be expressed as 

Li j2 (x), x eIq = (— oo, s{\; 

W t (x) = { ma x{V( Sl ),V(s i+1 )}, xel i = (s i ,s i+1 }; (IV) 

Lmt— l,mt 0*v> ^ ^ ^m t — (^rra t j ~)~C>o). 

The construction described above leads to the simplest proposal density possible, i.e., a collection of 
uniform pdfs with two exponential tails. Figure |4]^b) shows an example of the construction of the proposal 
using this approach. 

Note that, even when the target pdf is log-concave, using these two procedures the resulting ARMS- 
type techniques are not reduced to standard ARS, since there is no guarantee that the resulting proposal 
is above the target pdf. Indeed, when V{x) = log[p(a;)] is concave, the first procedure described 
by Eq. ( fTTj ) will produce a proposal which is always below the target except inside the first and 



last intervals. However, the second procedure described by Eq. ( fTTj ) can be easily modified to yield 
irt(x) = exp(Wt(x)) > p(x) = exp(V(x)) Vx G V for log-concave target densities, if the position of the 
mode is known or if an upper bound for V(x) is available and the location of the mode can be estimated 
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(a) (b) (c) 

Fig. 4. (a) Example of the construction of Wt (x) using the procedure described in Eq. (|16|. (b) Example of the construction 



of W t (x) using the procedure given by Eq. 1 17 1, (c) If the function V(x) is concave, the procedure in Eq. 1 17 i can be modified 



to obtain nt(x) > p(x) for all x £ T>. Indeed, this figure provides an example of the construction of Wt(x) when an upper 
bound B > V(x) is available and the location of the maximum can be estimated. 



(e.g. using the sign of the first derivative). Figure |4jc) depicts an example of the construction of Wt(x) 
when V(x) is concave using an upper bound B > V(x) to obtain TTt(x) > p{x) \/x £ V. Moreover, 
it is important to remark that the procedure in Eq. (F7\ , without any modifications, converges virtually 
to produce ir t (x) = exp(Wj(x)) > p{x) = exp(V(x)) almost everywhere when V(x) is concave and 
t — > +oo. 

Finally, we note that we could also apply the procedures proposed for adaptive triangular Metropolis 
sampling (ATRIMS) and adaptive trapezoid Metropolis sampling (ATRAMS) to build the proposal, even 
though the structure of these two algorithms is completely different to ARMS pCai et ah] 2008]. In both 



cases the proposal is constructed directly in the domain of the target pdf, p(x), rather than in the domain 
of the log-pdf, V(x) = log(p(x)). 

For instance, the basic idea proposed for ATRAMS is using straight lines, Lj j + i(x)J^] passing through 
the points (si,p(si)) and (si + i,p(si + i)) for i = 1, . . . , mt — 1 and two exponential pieces, Eq(x) and 
E mt {x), for the tails: 



7T t (x) OC 7Tf(x) = < 



E (x), x e J = (-oo,si]; 

L iti+ i(x), x eli = (si,s i+1 \, i = l,...,mt - 1; ( 18 ) 
E mt (x), x G I mt = (s mt , +oo). 



2 Note that we use the symbol ~ to distinguish Lij+i(x) and Li > i+i{x). Indeed, Li^+iix) is a straight line built directly in 
the domain of p(x) whereas the linear function Li,i+i(x) is constructed in the log-domain. 
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Unlike in | Cai et al. 2008], here the tails Eq{x) and E mt {x) do not necessarily have to be 



equivalent in the areas they enclose. Indeed, we may follow a much simpler approach calculating two 
secant lines £1,2(2;) and L mt ^ mt {x) passing through (s 1 ,V(s 1 )), (s 2 ,V(s 2 )), and (s mt -i, V(s mt -i)), 
(sm t ,V(s mt )) respectively, so that the two exponential tails are defined as Eq{x) = exp-fXi^ic)} and 
E mt (x) = exp{L mt _i )mt (x)}. 

Figure [5] depicts an example of the construction of TTt(x) using this last procedure. Note that drawing 



samples from these trapezoidal pdfs inside I{ = (s», Sj+i] can be easily done [Cai et al. 2008 Devroye 



1986[ . Indeed, given u', v' ~ U([s h s i+1 ]) and w' ~ U([0, 1]), then 

miniu',^'), w' < —, — P^r \\ 

r / n 1 \ p( s i) 

1 ' ' ' — P(Si)+p(Si +1 ) ' 

is distributed according to a trapezoidal density defined in the interval Zj 



(19) 



[Sj, Si+1 J. 




All these procedures reduce the computational cost of the so called ARMS-type algorithms. However, 
it is important to remark that they could be inefficient within the standard ARMS structure without the 



proposed enhancements, as shown in Section VII and due to the reasons explained in Section III An 
important exception is the simple procedure described by Eq. ( fP7| ), which leads to very good results 
also in combination with the standard ARMS (although even better results are obtained for A 2 RMS and 
IA 2 RMS), probably due to the fact that the proposal stays above the target almost everywhere. 
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A. Initial support points and tails 

Depending on the type of construction, the minimum possible number itiq of initial support points 



is at least 2 with the procedures described above or 3 with the procedure in [Gilks et al. 1995b |. In 
general, to avoid numerical problems the user must only assure that p(s') > for at least one support 
point s'. However, to speed up the convergence and render the algorithm more robust (avoiding any kind 
of numerical problems), a greater number of support points can be used. Moreover, observe that there is 
clearly a trade-off between the initial performance and overall computational cost. 

It is important to note that the dependence on the choice of the set Sq can be reduced adding more 
probability in the tails of the proposal PDF, initially. Indeed, consider the two slopes l t .\ and l t .2 of the 
linear functions that form Wt{x) in Xq and Z nit+ \, constructed as explained abovd^ For instance, we can 
add more probability in the tails of the proposal just by setting 

l' tA =l tA (l-f3e- at ), 

(20) 

4, 2 = Ml-/3e- a '), 

where < /3 < 1, a > 0, and /[ i — > l t i, i = 1,2, for t — > +oo. Moreover, the proposed algorithms work 
even if the target PDF has heavy tails on account of the second test in A 2 RMS and IA 2 RMS that allows 
the incorporation of new support points where p(x) > nt(x\St). However, if we know what kind of tails 
the target distribution has, we can speed up the convergence and reduce the burn-in period by changing 
the construction of the tails of the proposal. For instance, we could use a function of type log[l/a; 7 ], 
with 7 > 1, to build Wt(x) in Zo and X mf+ i, instead of straight lines. 

VII. Simulations 

In this section we compare the performance of the two proposed algorithms with the standard ARMS 
scheme using four alternative procedures to build the proposal pdf irt(x) in all cases: the one proposed by 
standard ARMS and described by Eq. ( [TT] ) (procedure 1), the two simpler procedures given by Eqs. ( fl6] ) 
and ( fTTj ) (procedures 2 and 3, respectively), and the procedure described in Eq. ( fT8] l, which is adapted 
from the ATRAM technique flCai et al.| |2008t (procedure 4)Q 



3 Note that, in general, a control on the slopes l t ,i and lt,i is needed to attain a proper proposal PDF since it is necessary 
that l t ,i > and l t ,i < 0. 

4 Matlab code using the procedure in Eq. {T7]l to build the proposal is provided at http : //a2rms . sourceforge . net/. 
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For the simulations, we consider a target pdf, p a (x) oc p(x), generated as a mixture of 3 Gaussian 
densities, 

p (x) = 0.3A/"(x; -5, 1) + 0.3A/"(x; 1, 1) + 0.4AA (x; 7, 1), (21) 

where N(x; fi, a 2 ) denotes a Gaussian pdf with mean \i and variance a 2 . In order to analyze the accuracy 
of the different algorithms, we estimate the mean of p {x) using the generated samples, comparing it to 
the true value, E{p Q (x)} = 1.6. Moreover, we also provide an estimation of the linear correlation among 



consecutive samples of the Markov chain, an approximation of the distance D n \ p (t) in Eq. ( pL0] > between 
the constructed proposal and the target pdf, the number of rejections and the number of support points 
at each iteration. 

We consider iV = 5000 iterations of the Markov chain and use an initial set Sq = {s\ = — 10, S2 = 
a, S3 = b, S4 = 10} formed by mo = 4 support points, where a and b (with a < b) are chosen randomly 
inside the interval [—10, 10] in all casesj^Moreover, in the first proposed scheme (A 2 RMS) we set K = N, 
i.e. we never stop the adaptation and use all the samples generated by the chain for the estimation. 



Tables III IV and [v] show the results obtained with the different algorithms (standard ARMS, A 2 RMS 



and IA 2 RMS), averaged over 2000 runs, using all the N = 5000 samples generated by the Markov 
chain and different procedures to build the proposal pdf. The columns of the Tables show (from left to 
right): the procedure used to build the proposal pdf, irt(x); the estimated mean (the true value is 1.6); the 
standard deviation of the estimation; the estimated linear correlation between consecutive samples of the 
chain; the average final number of support points]^] the average number of rejections in the first control 
(RS test); the average number of rejections in the second control (MH step) Q and an approximation of 
the discrepancy between the target and the proposal pdfs, given by D n \ p (t), after N = 5000 iterations. 

First of all, from these three tables we note that the two proposed algorithms provide better results 
than the standard ARMS in all cases, regardless of the scheme used to build the proposal. This can be 
seen by the slight improvement in the estimated mean averaged over the 2, 000 runs (we recall that the 
true mean is 1.6) and, most notably, by the large decrease in the standard deviation of the estimation. 

5 Note that, in general, it is necessary to check that the slope of the proposal is positive inside To and negative inside X mt in 
order to obtain a proper proposal pdf. Therefore, we fix si = — 10 and S4 = 10 just to guarantee this for t = 0, thus avoiding 
the problem of having to control the slopes of Wt(x) inside the first and last intervals, To = (—00, si] and T mt — [s mt , +00) 
respectively. Note also that support points smaller than -10 or larger than +10 can be added as part of the adaptive construction 
of the proposal without compromising the slope of the tails in any of the approaches. 

6 The number of support points is equal to the sum of the number of rejections in the two controls plus the 4 initial points. 



This test does not exist in the standard ARMS structure. Hence, the value of this column in Table 



III 



is always 0. 
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TABLE III 

Results with N = 5000 for the standard ARMS structure 



Procedure 


Estimated 


Std of 


Estimated 


Avg. num. 


Avg. num. of 


Avg. num. of 


Approx 


to build 


mean 


the 


correlation 


of support 


rejections in 


rejections in the 


of 


■ft(x) 




estimation 




points 


the RS test 


second control 


D n \ p (t) 


1 


1.6480 


0.7301 


0.3856 


65.8717 


61.8717 





3.0020 


2 


1.7541 


1.0908 


0.7722 


12.0492 


8.0492 





8.0516 


3 


1.5935 


0.2300 


0.6132 


164.1881 


160.1881 





6.1518 


4 


1.5670 


0.4961 


0.7083 


37.8231 


33.8231 





7.1339 



TABLE IV 

Results with N = 5000 for the A 2 RMS structure 



Procedure 


Estimated 


Std of 


Estimated 


Avg. num. 


Avg. num. of 


Avg. num. of 


Approx 


to build 


mean 


the 


correlation 


of support 


rejections in 


rejections in the 


of 


■ft(x) 




estimation 




points 


the RS test 


second control 




1 


1.6244 


0.1184 


0.0038 


94.4969 


80 


10.4969 


0.0613 


2 


1.7381 


0.2258 


0.0229 


86.1017 


9.8948 


72.2069 


0.0580 


3 


1.5984 


0.0797 


0.0026 


317.4166 


305.0933 


8.3232 


0.3110 


4 


1.6051 


0.0854 


0.0060 


92.3517 


54.8663 


33.4855 


0.0590 



TABLE V 

Results with N = 5000 for the IA 2 RMS structure 



Procedure 


Estimated 


Std of 


Estimated 


Avg. num. 


Avg. num. of 


Avg. num. of 


Approx 


to build 


mean 


the 


correlation 


of support 


rejections in 


rejections in the 


of 


■Kt{x) 




estimation 




points 


the RS test 


second control 


D n \p(t) 


1 


1.6233 


0.1238 


0.0041 


94.8435 


81.1630 


9.6805 


0.0609 


2 


1.7244 


0.2194 


0.0203 


85.6441 


10.0169 


71.6271 


0.0565 


3 


1.6007 


0.0950 


0.0021 


317.5360 


306.2528 


7.2832 


0.3009 


4 


1.6011 


0.1308 


0.0054 


92.1333 


55.0159 


33.1175 


0.0582 



This can also be clearly appreciated in Figure [6] where the average value of the estimated mean, plus 
and minus the average standard deviation, is plotted for t = 0, 1, . . . , 5000. Although all the curves 
converge to the true mean, the ones obtained using the two variants of ARMS proposed (A 2 RMS and 
IA 2 RMS) show clearly a reduced variance in the estimation, especially after the 5000 iterations have 
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been performed. 



Regarding the comparison of ARMS with the two proposed algorithms performed in Tables III IV 
and [Vj we also notice that the correlation after N = 5000 iterations is much lower for the two variants 
of ARMS proposed than for the standard ARMS. Furthermore, the convergence of nt(x) to p(x) also 
improves greatly, as evidenced by the value of D n \ p (t), which is more than one order of magnitude lower 
than for the standard ARMS after N = 5000 iterations, both for A 2 RMS and IA 2 RMS. 

We remark that all these improvements come at the expense of a larger number of support points 
(whose amount remains moderate anyway, not increasing without bound as t — > oo), and only at a 
slight increase in the computational cost, since both variants are very similar to ARMS, inserting just an 
additional step to incorporate support points inside regions where 7r t (x) < p(x). 

We also notice that the third procedure for building the proposal pdf, described by Eq. ( fTT] ), always 
provides the best results (regarding the average mean and standard deviation), even better than the first 
procedure with the standard ARMS, which is the strategy proposed in the original paper. Indeed, the 
third procedure provides a good approximation of all the modes of the target pdf more quickly than the 
other procedures, obtained at the cost of using a larger number of support points. This is owing to the 
procedure 3 having a greater "disposition" to propose and then incorporate support points in different 
regions of the domain, since it is formed by pieces of uniform pdfs. However, stepped lines are not the 
best approximation of a curve (in this case, the shape of p(x)), hence the distance D^ p (x) with the 
procedure 3 decreases more slowly than the rest. 

The second procedure attains a good approximation of the target p{x) (D n \ p (x) close to zero) 
always using the smallest number of support points. However, it provides the worst results in terms 
of estimation. Finally, the fourth procedure arguably provides the best trade-off among computational 
cost, approximation of the target density and accuracy of the estimation. 

VIII. Discussion 

The efficiency of MCMC methods depends on the similarity of the proposal and the target densities. 
In fact, the relationship between the speed of convergence of MCMC approaches and the discrepancy 
between the proposal and the target pdf has been remarked in different works (see e.g. [ Holden[ 1998]). 



Hence, in order to minimize this discrepancy and speed up the convergence, many adaptive strategies have 
been proposed in the literature [Andrieu and Moulines, 2006, An drieu and Thorns] 2008 1 Cai et al. [2008 



Gasemyr| [2003] |Griffin and Walker| |20TT] |Haario et al^ [200T] |Holden et aL| [2009 ] . Indeed, the automatic 



construction of a sequence of proposal pdfs that are both easy to draw from and good approximations 
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to the target density, is a very important issue (and a difficult task) to speed up MCMC algorithms. 

In order to draw from complicated univariate target densities, one popular approach for building 
the proposal adaptively is adaptive rejection Metropolis sampling (ARMS), which is an important 
generalization of standard adaptive rejection sampling (ARS). ARMS was introduced to allow drawing 
samples from univariate non-log-concave distributions through the addition of a Metropolis-Hastings step, 
and reduces to ARS for log-concave target pdfs. Furthermore, it is important to remark that ARMS is often 
used to handle multivariate distributions within MCMC methods that require sampling from univariate full 
conditional distributions, such as Gibbs sampling [Rob ert and Casella| 2004], the hit-and-run algorithm 
[Beli sle et al.| |1993| and adaptive direction sampling | |Gilks et al} [T994]. 

In this work, we have first highlighted that the approach followed by ARMS to build the proposal 
can be decoupled from the procedure used to incorporate points to the support set. Moreover, we have 
pointed out a structural limitation of ARMS (support points can never be added inside regions where 
TTt(x) < p(x)) that may prevent the target from converging to the proposal. 

In order to overcome this drawback, we have proposed two alternative ARMS-type adaptive schemes 
that guarantee both the convergence of the proposal to the target density and the convergence of the 
Markov chain to the invariant distribution, including the statistical information provided by every sample 
x' G T> such that nt(x') < p(x'). The first one (A 2 RMS) is a direct modification of the ARMS procedure 
that incorporates support points inside regions where the proposal is below the target with a decreasing 
probability, thus satisfying the diminishing adaptation property, one of the required conditions to assure the 
convergence of the Markov chain. The second one (IA 2 RMS) is an adaptive independent MH algorithm 
with the ability to learn from all previous samples except for the current state of the chain, thus also 
ensuring the convergence of the chain to the invariant density |Gasemyrj 2003 1 Holden et al. 2009}, 
[ Liang et al.[ 2010[ Chapter 8]. The underlying idea behind both of the proposed schemes is measuring 
the discrepancy between the proposal and target densities through the computation of ratios of the two 
pdfs. 

We have shown through a numerical example that the two novel approaches proposed provide an 
improvement in performance w.r.t. standard ARMS: more accurate estimated mean, less variance in 
the estimation, less correlation between consecutive samples and better convergence to the target pdf. 
Regarding this last issue, it is important to remark that the two ARMS-type procedures introduced 
guarantee that the discrepancy between the proposal and the target pdf decreases quickly, implying that 
the proposal becomes closer and closer to the target density. Note that this is not ensured in other adaptive 



MH algorithms |Haario et al. 200 1| Liang et aL| 2010 1, which simply optimize certain parameters of 
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a proposal pdf with a fixed analytical form, whereas here we improve the entire shape of the proposal 
density, ensuring the convergence to the target pdf. 

The price paid w.r.t. standard ARMS is a slight increase in computational cost, due to the single 
additional step to allow for the incorporation of support points inside regions where 7r t (x) < p(x), and a 
moderate increase in the number of support points, which remains bounded anyway. Moreover, the effort 
to code and implement the novel algorithms remains virtually unchanged, since the new control step in 
A 2 RMS and IA 2 RMS does not need additional evaluations of the proposal and the target pdfs. 

Finally, we have also noticed that these novel techniques allow us to reduce the complexity in the 
construction of the sequence of proposals w.r.t. standard ARMS, thus compensating the slight increase 
in computational cost. Consequently, we have also proposed different simpler procedures to build the 
sequence of proposal densities and tested them with the standard ARMS, A 2 RMS and IA 2 RMS. This 
is an important issue as evidenced, for instance, by the fact that the work of [Me yer et al.[ 2008 1 is 



dedicated exclusively to an alternative construction of the proposal pdf within the ARMS technique. 

One of introduced procedures, constructed using pieces of uniform pdfs plus two exponential-type pdfs 
for the tails, is particularly interesting due to its simplicity and good performance even with standard 
ARMS. Indeed, it could potentially be used to design efficient ARMS algorithms to draw directly from 



multidimensional target distributions. Moreover, the procedure based on the work of |Cai et alj |2008 1 
that uses a piecewise linear approximation of the target pdf and two exponential tails, provides the best 
trade-off among computational cost, approximation of the target density and accuracy of the estimation. 
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Fig. 6. Estimated mean as function of the iteration index t, averaged over 2000 runs, for the ARMS (Figures (a), (d), (g) and 
(j)), A 2 RMS (Figures (b), (e), (h) and (k)) and IA 2 RMS (Figures (c), (f), (i) and (1)). Figures (a), (b) and (c) correspond to 
■nt(x) obtained using Eqs. \\ \\-\\2\. Figures (d), (e) and (f) to 7Tt(x) given by 1 16 1. Figures (g), (h) and (i) to 7Tt(x) given by 
(17). Fi gures (j), (k) and (1) are obtained using the procedure described in Eq. {T8}. 



